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Abstract 

We  will  compare  speckle  imaging  reconstruction  results  for  several  speckle  imaging  approaches.  *  In  particular,  we 
will  compare  and  contrast  four  methods:  1)  Knox-Thompson,  using  a  hidden  phase-finder  in  the  object  spectrum 
phase  reconstruction;  2)  Knox-Thompson,  using  a  phasor-based  phase  reconstruction;  3)  Bispectrum,  using  only 
two  bispectmm  planes;  4)  Bispectrum,  using  four  bispectrum  planes.  In  each  application  of  the  four  approaches, 
we  first  calculate  the  modulus  of  the  object  spectrum  using  a  Wiener-Helstrom  filter  to  remove  the  speckle  transfer 
function.  The  methods  then  differ  only  in  their  object  spectrum  phase  reconstructions.  In  the  simulations,  we  will 
assume  that  the  only  aberrations  are  those  introduced  by  atmospheric  turbulence,  setting  the  ratio  of  the  telescope 
diameter,  D,  to  the  Fried  Parameter,  ro,  equal  to  ten.^  Additionally,  we  assume  that  the  focal-plane  detector  array  is 
photon-noise  limited,  the  illumination  is  narrow-band  (essentially  monochromatic)  and  the  atmosphere  is  static 
during  the  each  data  frame.  First,  we  will  implement  all  four  methods  on  a  simple  binary  star  object  at  low  photon- 
per-frame  light  levels.  Next,  we  will  apply  the  methods  to  complex  extended  objects. 


Introduction 

Speckle  imaging  techniques  have  been  evolving  since  the  fundamental  idea  was  presented  almost  forty  years  ago.* 
The  critical  insight  demonstrated  how  moments  of  the  Fourier  transform  of  an  ensemble  of  short  exposures  contain 
information  out  to  the  diffraction  limit.  Many  variations  on  the  theme  have  been  implemented,  but  in  all  cases,  an 
ensemble  of  short-exposure  images  is  collected  and  then  post-processed  to  restore  the  object.  In  this  presentation, 
we  will  compare  speckle  imaging  reconstruction  results  for  several  speckle  imaging  approaches.  Specifically,  we 
will  compare  and  contrast  four  methods:  1)  Knox-Thompson,  using  a  hidden  phase-finder  in  the  object  spectrum 


1 


Report  Documentation  Page 

Form  Approved 

0MB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  0MB  control  number. 

1.  REPORT  DATE 

SEP  2009 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2009  to  00-00-2009 

4.  TITLE  AND  SUBTITLE 

Comparing  Speckle  Imaging  Methods 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PEREORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory, Directed  Energy  Directorate,Kirtland 
AEB,NM,87185 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR’S  ACRONYM(S) 

II.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

2009  Advanced  Maui  Optical  and  Space  Surveillance  Technologies  Conference,  1-4  Sep,  Maui,  HI. 


14.  ABSTRACT 

We  will  compare  speckle  imaging  reconstruction  results  for  several  speckle  imaging  approaches.  1  In 
particular,  we  will  compare  and  contrast  four  methods:  1)  Knox-Thompson,  using  a  hidden  phase-finder  in 
the  object  spectrum  phase  reconstruction;  2,3  2)  Knox-Thompson,  using  a  phasor-hased  phase 
reconstruction;  3)  Bispectrum,  using  only  two  hispectrum  planes;  4  4)  Bispectrum,  using  four  hispectrum 
planes.  In  each  application  of  the  four  approaches  we  first  calculate  the  modulus  of  the  object  spectrum 
using  a  Wiener-Helstrom  filter  to  remove  the  speckle  transfer  function.  The  methods  then  differ  only  in 
their  object  spectrum  phase  reconstructions.  In  the  simulations,  we  will  assume  that  the  only  aberrations 
are  those  introduced  by  atmospheric  turbulence,  setting  the  ratio  of  the  telescope  diameter,  D,  to  the  Fried 
Parameter,  rO,  equal  to  ten.5  Additionally,  we  assume  that  the  focal-plane  detector  array  is  photon-noise 
limited,  the  illumination  is  narrow-band  (essentially  monochromatic)  and  the  atmosphere  is  static  during 
the  each  data  frame.  First,  we  will  implement  all  four  methods  on  a  simple  binary  star  object  at  low 
photonper-  frame  light  levels.  Next,  we  will  apply  the  methods  to  complex  extended  objects. 


15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 

OE  PAGES 

1 9a.  NAME  OE 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

10 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


phase  reconstruction;  2)  Knox-Thompson,  using  a  phasor-based  phase  reconstruction;  3)  Bispectmm,  using  only 
two  bispectrum  planes  and  a  phasor-based  phase  reconstruction;"'  4)  Bispectmm,  using  four  bispectmm  planes  and  a 
phasor-based  phase  reconstmction.  In  all  applications,  we  first  calculate  the  modulus  of  the  object  spectmm  using  a 
Wiener-Helstrom  filter  to  remove  the  speckle  transfer  function.  The  methods  then  differ  only  in  their  object 
spectmm  phase  reconstmctions. 

The  first  method  solves  two-dimensional  difference  equations  for  the  phase  using  the  method  described  in  Reference 
3.  There,  we  demonstrate  that  the  object  spectmm  phase  can  be  decomposed  into  a  regular,  single-valued  function 
determined  by  the  divergence  of  the  phase  gradient,  as  well  as  a  multi-valued  function  determined  by  the  circulation 
of  the  phase  gradient;  this  second  function  has  been  called  the  "hidden  phase."  The  standard  least-squares  solution 
to  the  two-dimensional  difference  equations  always  misses  this  hidden  phase.  Reference  3  develops  a  solution 
method  that  gives  both  the  regular  and  hidden  parts  of  the  object  spectmm  phase. 

The  next  three  methods  all  use  phasor-based  phase  reconstmction  algorithms.  Here,  in  each  case,  we  develop 
iterative  improvement  algorithms  that  rapidly  converge  to  the  least-squares-best,  two-dimensional  phasor  array  for 
each  object  spectmm. 

In  our  applications,  we  will  implement  all  four  methods  on  a  simple  binary  star  object  at  low  photon-per- frame  light 
levels.  Then,  we  will  apply  the  four  methods  to  several  complex  extended  objects,  once  again  varying  the  photons 
per  frame.  In  the  simulations,  we  will  assume  that  the  only  aberrations  are  those  introduced  by  atmospheric 
turbulence  setting  the  ratio  of  the  telescope  diameter,  D,  to  the  Fried  Parameter,  ro,  greater  than  or  equal  to  ten.^  In 
all  cases,  we  apply  each  of  the  methods  to  aberrated  image  ensembles  of  two-hundred  short-exposure  data  frames. 

In  addition,  we  assume  only  photon  noise  in  the  short  exposures  while  neglecting  other  noise  sources.  We  will 
present  numerous  results,  describing  the  strengths  and  weaknesses  of  each  of  the  four  methods  applied  to  both 
simple  and  extended  objects. 


Reconstructing  the  Object  Spectrum  Modulus 

We  label  the  arrays  of  short-exposure  data  as 


d{l;i,j).  (1) 

The  index  /  =  1,2,  . . .  R  labels  the  frame,  while  the  next  two  indices,  (ij),  label  the  pixel  in  the  array.  We  take  a  two- 
dimensional  Finite-Fourier-Transform  (FFT)  on  the  pixel  indices  as 

d{l;m,n)  =  FFT  {d{l',i,j)}.  (2) 
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The  array  indices,  (m,  «),  label  the  elements  of  the  two-dimensional  FFT.  The  zero  spatial  frequency  point  is  labeled 
{m,n  )  .  We  then  form  the  unbiased  average  ^ 


C{m,n)  = 


)-(f  (/;»!,«)(=©(»!,«)  •H{m,n) 


(3) 


in  which  |©(OT,n)|  is  the  modulus  of  the  object  spectrum  (FFT)  and  H(m,n)  is  the  speckle  transfer  function;  this  can 

be  calculated  with  a  reference  star  or  even  calculated  analytically.  A  standard  Wiener-Flelstrom  deconvolution  filter 
then  yields  an  estimate  of  the  modulus-squared  of  the  object  spectrum.  In  the  following  comparisons,  we  use  the 
same  modulus  for  each  of  the  four  speckle  methods;  therefore,  they  only  differ  in  the  phase  reconstruction  methods. 


Reconstructing  the  Object  Spectrum  Phase 

The  object  spectrum,  essentially  the  FFT  of  the  object,  is  given  by  the  modulus  times  the  phase  factor  as 


©(ot,  n)  =  @(m,  n)  ■  exp(;  W{m,  n))  =  ©(w,  n)  ■  P(m,  n)  , 


(4) 


in  which  W{m,n)  is  the  object  spectrum  phase  array  and  P(m,n)  is  defined  as  the  object  spectrum  phasor  array.  In 
what  follows,  we  will  develop  and  compare  methods  that  use  a  Knox-Thompson  or  Triple-Correlation  data  base  to 
generate  phase  estimates.  Results  for  the  phase  array,  W,  or  phasor  array,  P,  will  then  complete  our  estimate  of  the 
object  spectrum.  In  each  case,  prior  to  applying  an  inverse  FFT  to  the  reconstructed  object  spectrum,  we  applied  the 
diffraction-limited  optical  transfer  function  for  the  telescope;  this  avoids  signal-to-noise  problems  at  the  higher 
spatial  frequencies. 


In  the  two  Knox-Thompson  methods,  called  KTA  and  KTB,  we  start  by  forming  the  unbiased  cross-spectrums  from 
an  ensemble  of  R  short  exposures 

R 


C^{m,n)  =  —  X  {d{l',m,n)d*[l;m  +  l,n)-  d’'[l;m  -i- 1,«)| 

^  /=i 


C^, [m,n)  =  —  X  {d{l’,m,n)  d' +  l)  -  d' -t  l)}. 

^  /=i 


(5) 


in  which  we  have  picked  one  unit  offset  in  each  direction  in  the  FFT;  other  offsets  can  easily  be  analyzed,  but  the 
unit  offset  choice  gives  the  best  accuracy  on  the  high-spatial-frequency  phase.  The  ensemble  averaging  over  many 
realizations  of  the  turbulence-induced  phase  screens  makes  the  cross-spectrum  transfer  functions  real-valued.  ^ 
Therefore,  the  turbulence-induced  phase  errors  are  averaged  away  and  the  phase  of  the  cross-spectrum  can  be 
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directly  related  to  phase  differences  in  the  object  spectrum.^  The  remaining  phases  of  these  two  cross-spectrums  lead 
to  two  two-dimensional  difference  equations  for  the  phase  of  the  object  spectrum  as 

C^{m,n)  =  I  C^[m,n)  \  exp(  i  V^)  =  \C  ^{m,n^P[m,n)P*  [m  -i- 1,  n) 

(6) 

C^,  {m,n)  =  I  Cj,  (m, «)  I  exp(  i  V^,  )=  |C^  [m,  «)j  P{m,n)P*  (m,  n  +  \) 
in  which  the  cross-spectrum  phases,  and  are  related  to  the  object  spectrum  phase,  W,  as 

V^{m,n)=  w{m,n)-  w{m  +  1,«) 

(7) 

Vy{m,n)=  w{m,n)-  w{m,n  +  l). 

In  the  first  method,  KTA,  we  solve  these  difference  equations  directly  for  both  the  regular  and  hidden  parts  of  the 
object  spectrum  phase  array,  W{m,n).  This  method  is  detailed  in  reference  3 .  There,  we  find  the  regular  phase  by 
minimizing  the  error  function 

^  ^  X|  I  Wg{m,n)+  Wg{m  -i- 1,«)|^  -i- 1  V^,{m,n)-  Wg{m,n)+  Wg{m,n  -i- 1)|  |  (8) 


while  the  hidden  part  of  the  phase,  the  part  driven  by  helical  phase  around  branch  points  in  the  object  spectrum,  is 
given  by 


W„{m,n) 


J_ 

2n 


^arctan 

kj 


n-l  -  .5 
m-k -.5 


Vx{k,l  + 1)-  vAk,l)- V^k  + 1,/)+  V^kj)] 


(9) 


Figs,  (la)  and  (lb)  demonstrate  the  importance  of  the  hidden  phase  component,  with  the  left-hand  reconstruction 
using  only  the  regular  phase,  while  the  right-hand  reconstruction  uses  both  the  regular  and  hidden  phase 
components;  obviously  the  hidden  phase  component  is  critical. 

For  the  second  method,  KTB,  we  once  again  start  with  Knox-Thompson  data.  In  this  case,  we  solve  for  the  phasor 
array  by  minimizing  the  least-squares  error  function  given  as 


-1  ^  'J  -1 

E=  ^  \C^(m,nj^P{m,n)P  {m  +  \,ri)-C^  (m,n)  -I-  X  \Cy{m,n)p{m,n)P  (m,  n  +  \)-  Cy  (m,  n) 


(10) 


m,n 


m,n 


A  variation,  with  respect  to  P(k,l),  leads  to  a  nonlinear  system  of  equations.  These  can  then  be  iteratively  solved  to 
convergence,  yielding  a  least-squares-best  phasor  array,  P(m,n) .  This  procedure  is  similar  to  the  method  used  by 
Fludgin  on  phase  differences.  ^ 
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In  the  two  bispectrum,  or  Triple-Correlation  methods,  referred  to  as  TCA  and  TCB,  we  start  by  forming  several  sets 
of  the  short-exposure  ensemble  averages  defined  by 


p{m,  n;m',n')- — ^  \d{l',m,n)d{l',m’ ,n  ;m  +  m',n  +  n')- noise  terms^  , 

R  1=1 


(11) 


in  which  (»!,«)and  range  over  all  spatial  frequencies  in  the  object  spectrum,  while  {m’,n')  are  fixed  array  elements. 
Each  value  of  the  pair  (m',  n')  defines  a  bispectmm  plane. 


In  TCA,  we  choose  two  bispectmm  planes  defined  by  m'  -m  +  l  ,n'  -n  and  m'  =  m,n'  =  «-i-l,in  which  {m,n) 


specifies  the  zero  spatial  frequency  element  in  the  spectram.  This  choice  is  essentially  identical  to  the  unit  shift  used 
in  Knox-Thompson.  Once  again,  the  object  spectmm  phasor  array  is  taken  to  minimize  the  least-squares  phasor 
error  given  by 


n;m  +  l,n^P{m,n)p{m  +  l,n)P  ^{m  +  l,n)- p{m,n\m  +  1,77 J 

m,n 

+  X  n;m,n  -i-  l)^P{m,n)p{m,n  -i-l)P  ^{m,n +  l)- /]{m,n;m,n  -t-l) 
m,n 


(12) 


Here,  once  again,  the  variation,  with  respect  to  P(k,l),  leads  to  a  nonlinear  system  of  equations.  These  can  then  be 
iteratively  solved  to  convergence,  yielding  a  least-squares-best  phasor  array,  p(m,  n) . 


The  constraction  of  TCB  is  based  on  four  closely-spaced  bispectmm  planes  with  [m'  -in  +  l ,  n  -n) , 

{m'  =  m  ,  n'  =  n  +  =  m  +  2  ,  n'  =  n),  and  {m'  =  m  ,  n'  =  n  +2).  Here,  the  least-squares  phasor  error  function 

depends  on  four  bispectmms  as 


1  ^ 

E  =  X  -I- l,n^P(m,n)p(m  +  l,n)P  (m  +  l,n)- /3(m,n’,m  +  l,n) 

m,n 

+  .  -I-  . 

-I-  X  n',m,n  +  2jp{m,n)p{m,n  +  2)P  ^(m,n  +  2)- p{m,n',m,n  +  2) 
m,n 


(13) 


Once  again,  the  least-squares-best  values  for  the  phasor  array  result  from  an  iterative  procedure  that  minimizes  this 
error. 
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Reconstruction  Results 

In  the  simulations,  we  will  assume  that  the  only  aberrations  are  those  introduced  by  atmospheric  turbulence,  setting 
the  ratio  of  the  telescope  diameter,  D,  to  the  Fried  Parameter,  ro,  equal  to  ten.^  We  could  generalize  the  simulations 
to  include  static  and  dynamic  telescope  aberrations.  Additionally,  we  assume  that  the  focal-plane  detector  array  is 
photon-noise  limited,  so  that  we  neglect  dark-current,  as  well  as  read  noise.  Also,  we  only  process  the  data  on 
(64X64)  pixel  elements.  We  fix  the  final  F/#  of  the  telescope  so  that  the  diffraction-limited  point-spread-function  is 
approximately  two,  to  three  pixels  in  diameter.  In  all  cases,  we  used  the  same  speckle-derived  modulus,  while  each 
phase  or  phasor  reconstmctor  used  the  algorithms  described  for  KTA,  KTB,  TCA  and  TCB  in  the  last  section.  In 
each  case,  we  applied  the  diffraction-limited  optical  transfer  function  for  the  telescope  prior  to  applying  an  inverse 
FFT. 

First,  we  implemented  all  four  methods  on  a  simple  binary  star  object  at  low  photon-per-frame  light  levels.  For  each 
case,  we  used  an  ensemble  of  two  hundred  short-exposure  frames  with  one  hundred  photons  per  frame.  Figs.  (2a, 

2b)  show  a  single  frame  diffraction-limited  binary  image  with  photon  noise,  as  well  as  the  long-term  aberrated 
image.  The  results  for  all  four  speckle  methods  are  shown  in  Figs.  (2c-2f).  They  are  all  acceptable  reconstructions 
of  the  simple  binary  object. 

Next,  we  implemented  all  four  methods  on  a  simple  four-star  asterism  at  low  photon-per-frame  light  levels.  For 
each  case,  we  used  an  ensemble  of  two  hundred  short-exposure  frames  with  five  hundred  photons  per  frame.  Figs. 
(3a,  3b)  show  a  single  frame  diffraction-limited  asterism  image  with  photon  noise,  as  well  as  the  long-term  aberrated 
image.  The  results  for  all  four  speckle  methods  are  shown  in  Figs.  (3c-3f).  In  this  case,  once  again,  the  methods 
offered  acceptable  reconstructions. 

Finally,  we  reconstructed  a  more  complicated  and  extended  object  at  higher- light  levels.  We  used  an  ensemble  of 
two  hundred  short-exposure  frames,  but  increased  the  photons  per  frame  to  ten  thousand.  Figs.  (4a,  4b)  show  a 
single-frame,  diffraction-limited  image  with  photon  noise  as  well  as  the  long-term  aberrated  image.  The  results  for 
all  four  speckle  methods  are  shown  in  Figs.  (4c-4f).  For  this  extended  object,  the  differences  between  the  four 
speckle  reconstruction  approaches  are  more  pronounced,  with  KTB  and  the  closely  related  TCA  giving  the  best 
results  and  TCB  the  worst. 


Discussion  and  Conclusion 

The  examples  demonstrated  reconstructions  on  200  short-exposure  frames  in  the  presence  of  significant  turbulence. 
For  simple  objects  with  object  spectrums  extending  out  to  appreciable  spatial  frequencies,  all  of  the  methods,  KTA, 
KTB,  TCA,  and  TCB  can  work.  For  extended  objects,  with  much  more  compact  object  spectrums.  Triple¬ 
correlation  methods  such  as  TCB  are  suspect.  ’  Methods  that  incorporate  higher- spatial  frequency  data  by 
incorporating  a  large  number  of  bispectrum  planes  will  encounter  difficulties,  since  the  extended  objects  typically 
have  a  small  object  spectrum  modulus  at  the  larger  spatial  frequencies.  When  we  recall  that  our  TCB  method  only 
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used  four  closely-spaced  bispectrum  planes  and  is  already  manifesting  problems  on  an  extended  object,  we  can  well 
appreciate  how  additional  bispectrum  planes  can  engender  even  more  reconstruction  problems. 

Although  our  Knox-Thompson  methods,  KTA  and  KTB,  as  well  as  the  closely  related  TCA,  worked  in  these 
examples  and  in  many  other  test  problems,  several  questions  remain:  How  will  the  methods  compare  when  more 
realistic  noise  sources  and  additional  telescope  aberrations  are  incorporated  into  the  simulations?  How  will  the 
methods  work  on  real  data?  How  can  we  know  that  our  phasor-based  error  functions,  Eqs.  (10),  (12)  and  (13),  are 
the  best  choices?  In  general,  is  there  a  truly  optimum  speckle  reconstmction  method,  an  optimum  method  for  both 
simple  and  extended  object  reconstructions  at  low-light  levels  and  significant  atmospheric  aberrations?  These 
questions,  along  with  other  evolving  issues,  should  be  addressed.  Based  on  the  present  set  of  reconstruction  results, 
the  phasor-based  Knox-Thompson  method,  KTB,  would  appear  to  be  close  to  the  optimum  method,  working  well  on 
both  simple  and  extended  objects. 
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Fig.  la.  Reconstruction  without  hidden  phase.  Fig  lb.  Reconstruction  with  hidden  phase. 
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Fig.  2a.  Binary  star,  diffraction  limited  image. 


Fig.  2b.  Binary  star,  long  term  aberrated  image. 


Fig.  2e.  Triple-Correlation  method  A.  Fig.  2f  Triple-Correlation  method  B. 


Fig.  3a.  Asterism,  diffraction  limited  image. 


Fig.  3c.  Knox-Thompson  method  A. 


Fig.  3e.  Triple-Correlation  method  A. 


Fig.  3b.  Asterism,  long  term  aberrated  image. 


Fig.  3d.  Knox-Thompson  method  B. 


Fig.  3f  Triple-Correlation  method  B. 


Fig.  4d.  Knox-Thompson  method  B. 
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Fig.  4f.  Triple-Correlation  method  B. 
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